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We present a theoretical study of the ground state of a Bose-Einstein condensate with repulsive 
inter-particle interactions in a double-well potential, using a restricted variational principle. Within 
such an approach, there is a transition from a single condensate to a fragmented condensate as 
the strength of the central barrier of the potential is increased. We determine the nature of this 
transition through approximate analytic as well as numerical solutions of our model, in the regime 
where the inter-particle interactions can be treated perturbatively. The degree of fragmentation 
of the condensate is characterized by the degrees of first-order and second-order spatial coherence 
across the barrier. 
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I. INTRODUCTION 

The recent experimental demonstration of interference 
phenomena in Bose-condensed atomic gases [Q, moti- 
vates a study of the spatial coherence of a condensate 
in a double-well potential. In particular, we are inter- 
ested in the loss of spatial coherence that can occur at 
zero temperature due to fragmentation of the conden- 
sate. A fragmented condensate is one for which there 
is a macroscopic occupation of two or more orthogo- 
nal single-particle wavefunctions. If the occupied single- 
particle wavefunctions are spatially well separated, co- 
herence over the spatial extent of the entire system will 
be lost, persisting only over the spatial extent of each 
fragment. 

As Nozieres has pointed out, for repulsive inter- 
particle interactions it is the exchange energy that typ- 
ically prevents fragmentation into a number of degener- 
ate (or nearly degenerate) single-particle wavefunctions. 
However, this argument is inapplicable for bosons in an 
external potential with several local minima, since single- 
particle wavefunctions that are localized about these 
minima may have very little overlap with one another, 
thereby leading to a very small exchange energy. More- 
over, since the self- interaction energy in such a frag- 
mented condensate is smaller than that of a single con- 
densate, it is possible for the total interaction energy to 
be smaller as well. Although every particle in the frag- 
mented condensate will pay a price in kinetic energy to 
occupy localized wavefunctions, the overall energy may 
still be less than that of a single condensate. Indeed, it 
can be shown that in the limit of a symmetric double- 
well potential with an infinitely strong central barrier, 
one can always find a fragmented state that has a total 
energy lower than any single condensate ||]. 

Thus we have the following situation in a double- 
well potential: in the absence of any central barrier the 
ground state is well approximated by a single conden- 
sate, while in the presence of an infinitely strong barrier 



it is well approximated by a fragmented condensate. It 
is clear therefore that there must be a transition between 
these two extremes as one increases the strength of the 
barrier. The first goal of this paper is to propose a theo- 
retical model for describing this transition. Specifically, 
we argue for an approximation of the fully interacting 
ground state that is more general than a Fock state, and 
that can be said to describe 'partial fragmentation' of the 
condensate. The equations that such a state must satisfy 
are derived within a variational approach. The second 
goal of the paper is to solve these equations in a regime 
where the inter-particle interactions can be treated per- 
turbatively. Numerical solutions of the equations and 
analytic approximations to these solutions are obtained 
within this regime. It should be noted that this limit 
is inappropriate for the description of the MIT conden- 
sate interference experiment Q|, and consequently, our 
results do not specify the nature of the transition for this 
experimental set-up. Nonetheless, we expect the generic 
features of the transition to persist in the experimentally 
relevant regime. 

We pause to consider previous treatments of this topic 
and their relation to this work. Rohrl et al. Q have pro- 
vided a model of the MIT condensate interference exper- 
iment; however, it is a mean field analysis and therefore 
cannot describe fragmentation. Fragmentation in the 
case of attractive inter-particle interactions has been con- 
sidered by Wilkin, Gunn and Smith but this effect is 
qualitatively different from that of the repulsive case. Fi- 
nally, Milburn et al. Q have considered the energy eigen- 
states of a Bose-Einstein condensate in a double-well po- 
tential, and predict fragmentation when the inter-particle 
interactions are sufficiently strong. However, these au- 
thors consider only traps with weakly coupled wells, and 
therefore cannot determine the degree of fragmentation 
of the ground state in the regime of low barrier strengths 
where the coupling between the wells is strong. More- 
over, this paper does not address the issue of the spatial 
coherence of the ground state. Since the presence of long- 
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range order is a defining characteristic of Bose-Einstein 
condensation, it is critical to understand the manners in 
which this spatial coherence can be lost. 

There has also been theoretical work on the problem 
of Bose condensates containing atoms in two different in- 
ternal states, which is analogous to the two well problem. 
Esry et al. || determine the probability distributions and 
lifetimes of two interacting condensates in different inter- 
nal states confined to the same trap, while Cirac et al. ||] 
as well as Steel and Collett || consider the ground states 
of such a system when the internal atomic states can 
be controlled by a Josephson-like laser coupling. In this 
case, the distinguishability of atoms in the two conden- 
sates is ensured by their internal state rather than their 
single-particle wavef unctions. However, the requirement 
of orthogonality in the spatial degrees of freedom is indis- 
pensable in a multiple well problem since it is precisely 
the shape of the single-particle wavefunctions that deter- 
mines the degree of fragmentation in the system. 

The remainder of the paper is organized as follows. In 
section II we present our model and define some useful 
measures of spatial coherence of the condensate. In sec- 
tion III we present approximate analytic solutions of our 
model in the regime of nearly non-interacting particles, 
and compare these to a numerical solution for a partic- 
ular choice of the external potential. The experimental 
signature of fragmentation and finite temperature effects 
are discussed in section V, followed by our concluding 
remarks in section VI. 



II. THE MODEL 

A. The basic approach 

Our system consists of an even number, N, of spin- 
less bosons at zero temperature. We model the inter- 
actions by a two-particle pseudopotential in the shape- 
independent approximation, V^r, r ) = g5(r — r ) with an 
interaction strength g = Aira sc h 2 /to, where a sc is the s- 
wave scattering length, and to is the mass of the bosons. 
The Hamiltonian is given by |10| 



H 



d 3 i 



-|^t(r)V 2 #(r) + t/(r)#t( r )i|r( r ) 
+ §# t (r)*t(r)*(r)*(r) 



(1) 



where ^(r) is the quantum field operator, and U(r) is 
the external potential. The external potential is taken to 
exhibit a single minimum along the y and z axes, and a 
double minimum along the x axis. It is also taken to be 
symmetric about x = 0. 

In order to capture the phenomenon of fragmentation 
in our model of the ground state, we must go beyond a 
mean field analysis. Specifically, we consider arbitrary 



superpositions of Fock states where up to two single- 
particle states are occupied. This corresponds to pos- 
tulating a state vector of the form 



JV 

£ 

ATi=0 



C Nl |iV X) iV 2 ) WiA)) 



(2) 



where 



\N U N 2 



(a\)* (4)^ 



\vac) 



is the Fock state in which Ni particles occupy the single- 
particle state <pi (r) = (r| a 



(r ak \vac) 



x \vac) and N 2 particles oc- 
The total number of particles 



cupy (fc (r) 

is fixed, N2 = N — N\, the vector consisting of the set 
of coefficients Cjvj is normalized, and the single-particle 
wavefunctions 0i an d 02 are both normalized and or- 
thogonal to one another. 

The state (||) is certainly not the most general state 
one can consider. Indeed, a state such as (||) would be a 
poor choice if one were interested in studying the deple- 
tion of a single condensate due to interactions, since one 
there expects a certain fraction of the particles to be dis- 
tributed among a macroscopic number of single-particle 
states. However, in this paper we are interested in the 
possibility of the particles being redistributed into a few 
single-particle states that are each macroscopically occu- 
pied. We restrict ourselves to two single-particle states 
because the double- well geometry we are considering en- 
courages fragmentation into two pieces ||. Although 
it may be energetically favorable to fragment into more 
than two pieces at very high particle densities, we defer 
consideration of this possibility to a later work. 

Among the many-body states defined by (||), we con- 
sider only those which have the same symmetry as the 
Hamiltonian under reflections about x = 0. This im- 
plies that the single-particle wavefunctions are mirror im- 
ages of one another across x = within a phase factor, 
4>i(— x, y, z) = e <fi 2 (x, y, z), and that the coefficients sat- 
isfy CjVi = Cjv-jVi- With this assumption, and choosing 
0i and 02 to be real, the Hamiltonian takes the form 

H 2 = e n N + (e 12 + gT^N - 1)) (a\a 2 + a\a^j 
(jvf + NI-N 
+^-^ L {a\a\a 2 a2 + a\a\aiai + 4NiN 2 ), (3) 



where iVi = a[a\,N 2 = a 2 a 2 , N = N\ + N 2 , and where 
e n = J d 3 rMr) (-^V 2 + C/(r)) 1 (r), 
ei2 = J d 3 rMr) (~|^V 2 + U(r)\ 2 (r), 



To 



^0i(r), 
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Ti= J d 3 r^(r)0 2 (r), 
T 2 = / d 3 r^(r)02( r ). 



These quantities have the following physical interpreta- 
tion: exi is the single-particle energy for the state <j>x; t\ 2 
is proportional to the inversion frequency of a single par- 
ticle in the external potential; finally, To quantifies the 
self-interaction energy, while T± and T 2 both quantify the 
cross-interaction energy. 

In order to facilitate comparison of our work with ear- 
lier studies ||, we re-express the Hamiltonian in 
terms of operators satisfying angular momentum com- 
mutation relations, rather than in terms of the creation 
and annihilation operators we have employed thus far. 
We introduce the operators 



{a\a\ — a\a 2 ), 



which form an angular momentum algebra with total an- 
gular momentum j = N/2 [Q]. In terms of these opera- 
tors, the Hamiltonian can be rewritten as 



H 2 = E + 2 (e 12 + gTi(N — 1)) J z + 2gT 2 J 2 z . 

'x ' 



-g(T -T 2 )4 



where 



E = e n N + In{N + 2){gT + gT 2 ), 



and where we have used J 2 + Jy+J z — j(j + 1) to elimi- 
nate Jy from the expression. The observable correspond- 
ing to J x is the particle number difference between the 
localized states 4>\ and <j) 2 . Defining the wavefunctions 
cj> s = 2~ 1 /2(0 1 + 2 ) and o = 2- 1 /2(^ 1 _ 2 ) ; which are 

respectively symmetric and antisymmetric about x = 0, 
one sees that J z can be rewritten as \{a\a s — a\a a ), 
where a\ and a\ are the creation operators associated 
with (j> s and <p a respectively. Thus, J z corresponds to 
the particle number difference between the symmetrized 
states 4> s and <fi a . Finally, J y corresponds to the conden- 
sate momentum. Since we are considering the ground 
state it follows that (J y ) — 0, and since the ground state 
is symmetric under reflections about x = it follows that 
( J x ) = 0. 



B. A restricted variational principle 

We now turn to the problem of identifying the ground 
state of our model. This can be achieved by minimizing 



the expectation value of H 2 with respect to variations in 
both the coefficients Cjvi and the single-particle wave- 
function <pi, subject to the constraints that the set of 
coefficients is normalized, and §\ is both normalized and 
orthogonal to <f> 2 , its mirror image about x — 0. However, 
it is in fact more convenient to minimize the expectation 
value of H 2 with respect to variations in <p s and 4> a rather 
than 4>i- The reason for this is that no constraint corre- 
sponding to orthogonality is required when working with 
cj) s and <f> a , since they are orthogonal by construction; as 
a result the analysis is simplified. 

We begin with the variation of <p s and <f> a , implement- 
ing the normalization constraints through Lagrange mul- 
tipliers E s and E a respectively. This results in two cou- 
pled non-linear Schrodinger equations for <p s and <j) a 



^-f/(r)+ 5 r°^ (r) 
-<?r^(r) 



Ira 



<f> a (r) =E a <j> a (r), (5) 



where 



T °a = (i a l a a) 2 - aia a ) j (ala a ) , and (6) 
r S = («l a Q a /3 a /3 + a^a^a a a a + 4aJ l a Q a| 3 a / 5^> / (a) a a a ) , 

and where the indices (a, /?) take the values (s, a) and 
(a,s). 

We now minimize the expectation value of H 2 with 
respect to variations in the Cn x , and implement the nor- 
malization constraint on the Cn 1 through a Lagrange 
multiplier E. This results in a five-term recurrence rela- 
tion for the coefficients 



Neu 



gT 



N) + 2gT 2 N 1 N 2 - E 



+ [e 12 + gT 1 (N-l)] 



y/ NM + ljC ^-! 

+ ^/N 2 (N 1 + l)C Nl+1 



gT 2 
2 



y/ fNj - l)JVi(JV 2 + l)(iV2 + 2)C Nl - 2 

+ ^{N 2 - l)JV 2 (iVi + l)(JVi + 2)C Nl +2 



(7) 



= 0, 



for each value of N\ ; E is immediately identified as the ex- 
pectation value of H 2 . The latter set of equations forms 
a matrix eigenvalue equation for the TV-element vector of 
coefficients C^ 1 . Given values for en, 612, Tq, T\ and T 2 , 
we can solve this equation by diagonalizing an N x N ma- 
trix with non-zero entries along five diagonals, a problem 
which is numerically tractable if the number of non-zero 
coefficients is not too large. 

Since equations (||) and (0) form a coupled set of equa- 
tions for <fi s , (j> a and the Cn x , we must in general solve 
these self-consistently. Of the many solutions thus ob- 
tained, the ground state is the one which minimizes the 
value of E. However, it is not obvious that the solution 
that minimizes E s and E a also minimizes E; thus it may 
be necessary to compare the energies of many solutions 
in order to find the ground state. 
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C. The regime of nearly non-interacting particles 

The full problem outlined above is rather complex. In 
this paper, we consider only perturbative solutions of Eq. 
(|J) in the nearly non-interacting regime, which we here 
define as the regime where the interaction energy is small 
compared to the difference between e« , the energy of the 
first symmetric excited state of the external potential, 
and e s ,the energy of the ground state. This is ensured by 
the criterion 



gNT « e 



(1) 



(8) 



Since To is on the order of the inverse of the volume of the 
trap, this criterion places an upper limit on the density 
of the condensate. 

In this regime, we can treat the non-linear terms in (|^) 
perturbatively. To obtain the expectation value of H 2 to 
first order in the perturbation, we need only solve for the 
eigenfunctions of (|5|) to zeroth order. Thus, we need only 
solve the two linear Schrodinger equations 



ft 2 V 2 
2m 



U(t) 



K (r) =0, 



(9) 



for a = s, a. In this case, the wavefunctions 4> s and <f) a 
are simply the two lowest single-particle energy eigen- 
functions of the external potential, and the assumption 
of a state of the form of (||) corresponds to a two mode 
approximation. 

The solutions of (^) determine the magnitudes of 
en, £12, To, T\ and T2, and these subsequently define the 
form of the recurrence relation ([?]) that must be solved 
to obtain the coefficients Cm x - Although the full results 
are presented in section III, it is illustrative to consider 
the ground state of H 2 in two particularly simple limits: 
that of no barrier and that of an infinitely strong barrier. 

In the absence of any barrier, |ei2 j ~ ei — e s , and since 
\Ti\ and T2 are on the order of To or less, it follows from 
the criterion (§) that |e 12 1 > NgT , Ng |T X | , NgT 2 . If 
we provisionally assume that the ground state fulfils the 



conditions that N 



T 



> J. 



Jl > , then we are led 



to approximate the Hamiltonian by 
H 2 — E + 2ei 2 J z . 



(10) 



The ground state of ( [l0| ) is simply the Fock state \N), 
which describes N particles occupying the single-particle 
ground state <p s . Since this solution satisfies our provi- 
sional assumption, the approximation is consistent. Such 
a Fock state is of course what one would expect for the 
ground state of a single well in the limit of nearly non- 
interacting particles. When this state is written in the 
form of (||, that is, in the basis of \N%, Nz),^ ^ states, 
rather than the basis of \N s ,N a )^ ^ ^ states, the coef- 
ficients C/Vj form a binomial distribution over N\, cen- 
tered at N/2. It seems appropriate to refer to any state 



of the form \N)^ o for macroscopic N and arbitrary </>o, 
as a 'single condensate'. In this paper, we are concerned 
only with single condensates wherein the single-particle 
wavefunction <pQ is symmetric about x = 0. We do not 
introduce any additional terminology to distinguish such 
a state from one with arbitrary </>o, since no confusion is 
likely to arise. 

In the limit of infinite barrier strength, the amplitudes 
of (f> s and 4>a at x = are necessarily zero, while <ps and 
4> a at x 5^ satisfy the same equation. Consequently, 
4> s and (j) a differ only in their symmetry under reflection 
about x = 0, and £12 = T\ = T 2 = 0. The Hamiltonian 
of (0) then reduces to 



H 2 =E + gT J 2 x 



(11) 



The ground state is \N/2,N/2}^ i ^ » , which describes 
two independent condensates, or in other words, a con- 
densate which is fragmented into two pieces. Since we 
are considering a potential well that is symmetric about 
x = 0, the two fragments are equally populated. It 
seems appropriate to refer to any state of the form 
\Ni 7 N 2 )( l p i ^ where Ni and N 2 are macroscopic and 
tfii and cf> 2 are orthogonal, as a 'dual condensate' pj| . In 
this paper we will be concerned only with dual conden- 
sates wherein N± — N 2 = N/2, and tf>\ and <p 2 are mirror 
images of one another across x = 0. 

The analysis above confirms, for the limit of nearly 
non-interacting particlws, the results of an earlier study 
Q: the ground state is well approximated by a single 
condensate at zero barrier strength, and a dual conden- 
sate at infinite barrier strength. At intermediate barrier 
strengths, we keep all the terms in H 2 for our calcula- 
tions. Although the cross-interaction terms are typically 
found to be small for generic shapes of the double-well 
potential, it is not obvious that these terms are negligible 
for an arbitrary potential, and thus we include them in 
our analytic results wherever possible. 



D. Measures of the degree of fragmentation 

Finally, in order to facilitate the interpretation of our 
results we highlight some observables that reveal the de- 
gree of spatial fragmentation of the condensate. The 
most useful observables for this purpose are those that 
probe the spatial coherence of the condensate across the 
barrier. In analogy to measures of optical coherence 
|p^| , we normalize the first-order correlation function, 

/3i(r,r') = /^(^"^(r'jV to obtain the degree of first- 
order spatial coherence between points r and r', 



5 «(r,r') = 



Pi(r,r') 



[pi(r, r)pi(r',r')] ! 



(12) 



Considering points r = (x, y, z) and r' = (—x, y, z) where 
x is positive and chosen to be sufficiently large so that 
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|0i (r)| < \<p 2 (r)\ and |0i(r')| > |0 2 (r')| , for any state 
of the form of g that is symmetric under reflection about 
x = 0, the quantity g^^(r,r') is in fact independent of r 
and r', and has the value 



N 



(13) 



We refer to simply as the degree of first-order spa- 
tial coherence across the barrier. It is straightforward to 
verify that it attains its maximum value of 1 for a single 
condensate and a value of for a dual condensate. 
The second-order correlation function p 2 (r, r') = 

^ t (r)* t (r')#(r')#(r)^ , which is simply the normally 

ordered density-density correlation, can be normalized 
to obtain the degree of second-order spatial coherence 
between points r and r', 

p 2 (r,r') 



,(2) 



(r,r') = 



[p 2 (r,r)p 2 (r',r')} 1 



(14) 



Defining C*- 2 -* in a manner completely analogous to 
we find 



C< 2 > = 



4(^) 2 



N f1 l Jf i 



(15) 



where AAq = ((^Nfj — ) 1 ^ 2 is the variance in the 

number of particles occupying the localized state 4>i . We 
refer to C^ 2 ' as the degree of second-order spatial coher- 
ence across the barrier. The variance in iVi for a single 
condensate {N}^ is that of a binomial distribution over 

Ni, namely \fN /2. For a dual condensate the number of 
particles in a well is fixed, so that AN\ = 0. As a conse- 
quence. C( 2 ) - 1 = for a single condensate \N) ^ , and 

— 1 = 2/N for a dual condensate. Since ANi is suffi- 
cient to specify C^ 2 ), while being simpler to interpret, we 
use AiVx together with to characterize our results. 



III. ANALYTIC APPROXIMATIONS AND 
NUMERICAL SOLUTIONS 

For a given shape of the double-well potential, it is 
straightforward to obtain the single particle ground state 
and first excited state by solving the linear Schrodinger 
equation (||). The localized single-particle states 4>\ and 
4>2 are simply the sum and difference of the single-particle 
ground and first excited states. Using these wavefunc- 
tions, the coefficients Cm x can be obtained by solving the 
recurrence relation (0). 

We here present approximate analytic solutions for the 
CjVi given an arbitrary shape of the double-well poten- 
tial. Subsequently, we consider a particular form of the 
external potential for which the single particle ground 
and first excited states are obtained numerically. This 
allows us to obtain a numerical solution for the C n x and 
to compare this solution to the analytic approximations. 



A. Continuum approximation 

Suppose the coefficients for the ground state satisfy the 
condition 

|Cjvi+i — CjvJ "C CjVi- 

It is then useful to construct a function C(u), defined 
over the real numbers, such that C(u) = Cn 1 at the 
discrete points u = N~ 1 {N\ — N/2), and such that 
C'(u) <C NC(u). Given this assumption of smoothness, a 
coefficient of the form Cn 1 + p , where p is a small integer, is 
well approximated by a Taylor expansion of C{u + p/N) 
to second order in p/N . In this way, the recurrence 
relation ((jj) for the Cn x can be recast as a second order 
differential equation for the function C(u). Moreover, any 
sum over Ni can be approximated by an integral over it. 
In particular, the constraint of normalization for the co- 
efficients is replaced by the constraint that the integral of 
C 2 (u) over all u is 1/N. The assumption of smoothness 
is readily verified to be appropriate for a single conden- 
sate, and we therefore expect it to continue to hold for 
solutions over a range of small barrier heights. 

We also make use of the fact that the coefficients Cn ± 
are significant only in the region where Ni < \/N, or 
equivalently, that the function C{u) is significant only 
where u < y/T/N. This follows from the fact that any 
set of coefficients that has significant amplitude outside 
the range N\ < \fN also has an energy that is larger 
than a single condensate; the interaction energy is greater 
since it scales with AN\, and the single-particle energy 
is greater since it is a minimum for a single conden- 
sate. It is therefore appropriate to expand each of the 
-/Vi-dependent terms as a power series in u : 



N 



x /N 1 (N 2 + l) = J2l n u n 1 



j^y/(N 1 -l)N 1 (N 2 + l)(N 2 + 2) = J nu n , 



which implicitly defines the /„ and J n . 

The second-order differential equation for C(u) we 
obtain is found to have a first-order term which 
can be eliminated by the substitution C{u) — 
C{u) exp(— a 2 nU 2n ) with an appropriate choice of the 
constants a 2n . It then follows that C(u) satisfies a second- 
order differential equation identical to that of a particle 
with position a in a one-dimensional potential well of 
even powers of u. Since the solution is only significant 
in the range u < y/l/N, we make the approximation 
that terms in the potential that are quartic or of higher 
order in u are small perturbations upon the quadratic 
term, and can be neglected. In this case, the function 
C(u) = C(u) exp(— tu 2 /2) satisfies the equation for the 
modes of a simple harmonic oscillator: 

- C"(u) + [y 2 + t 2 )u 2 C{u) = (rj + t)C{u), (16) 
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where, after expanding the I n and J n to leading order in 
l/N, 



r\ = 2N- 



,E/N - [e n + e 12 + Ng(±T + T x + |T a )] 



-ei2- Ng{T 1 +T 2 ) 



2N 2 



-ei2- Ng(T 1 +T 2 ) 

2 (1 - $0 (-ei2 - iVffTi) 
+iV ff (T + (4^ - 3) T 2 ) 



, and 



r = 2. 



(17) 



Thus, the solution for C(u) that minimizes 77, thereby 
minimizing the energy E, is a Gaussian 



G(«) = 



1 



1 



N (2tt 



with width 



2VN y/-A + y/A?TB' 



(18) 



(19) 



where 



A = 



B 



3 -e^-Ng^ 



2N - ei2 -JVfl(Ti 

-612 



T 2 ) 



, and 



- jVg(Ti + §T 2 - ±T ) 
-ei2 - JV ff (Ti + T 2 ) 



and with energy 



E = N 



47v^2 h (,£11 + £12) 

+ iV 5 (ATo + T 1 + fT 2 ) 



(20) 



The variance in N% for this solution is simply Na, while 
the degree of first-order spatial coherence is given by 



"(1 



-2a 2 ). 

absence 
NgT 2 , 



(21) 



In the absence of 

any barrier, |eia| 3> NgTo, Ng \Ti\ , NgT 2 , and it can be 
verified that eq. (19) predicts the appropriate value for 
this limit, namely ANi ~ l/2yN, the value for a sin- 
gle condensate. As the barrier strength is increased, the 
magnitudes of ei 2 ,Ti and T 2 decrease while the magni- 
tude of To does not vary significantly; it therefore follows 
from ([l9]) that AN\ will decrease with barrier strength. 
When AN\ falls below 1, the assumption of smoothness 
breaks down. Thus the range of validity of the continuum 
approximation is ANi > 1, or equivalently, C' 1 ' > 0.88. 

If the potential is such that Ng \T\ | , NgT 2 <C | e i2 1 con- 
tinues to hold as the barrier strength is raised from zero, 
then to leading order in l/N the expression for a simpli- 
fies to 

1 1 



2VN 



9T0 
2 F77i7 



depending only on the ratio of the interaction energy to 
the splitting between the symmetric and antisymmetric 
levels of the trap. 



B. Two-coefficient approximation 

At large barrier strengths, we make use of the following 
conditions: 



I7I <C 1, where 



f (f + l)(- £l2 - g^N-m) 



and 



|C| <C 1, where 



9T0 



N 2 gT 2 



(22) 



f (f + l)(-e 12 - g^N-m) 



(23) 



The first of these conditions is always satisfied for suffi- 
ciently strong barriers since in the limit of an infinitely 
strong barrier, ei 2 = T\ = 0, while To is finite. More- 
over, we have numerically verified that the second con- 
dition holds at sufficiently large barrier strengths for a 
variety of external potentials. Within the domain of ap- 
plicability of these conditions, we seek coefficients 
that satisfy the recurrence relation (^) to first order in 7 
and £. Dividing the recurrence relation by gTo, one finds 
that all terms involving N 2 gT 2 have a magnitude on the 
order of 7C and can therefore be neglected. In this limit, 
the following set of coefficients are a solution: 



C Nl = y/l - 2 7 2 for Ni = f 
= 7 for Ni = f 

= otherwise. 



1 K - l 



(24) 



I A (N — 2) _ ^ 2 ^ (25) 



The energy in this case is 
E = Ne n + I 



We refer to this approximation as the two- coefficient ap- 
proximation, and we dub any state of the form of ( j24| ) 
a perturbed dual condensate. For such a state, C' 1 ' and 
AN± are given by 



C {1) =27 V /l- 2 7 2 A /l+' 
V27. 



j?, and 



AN 



(26) 
(27) 



Keeping terms to first order in 7, and to leading order in 
powers of 1 /N, we haveC^ = 2 7 . The range of validity of 
the two-coefficient approximation is the range of barrier 
strengths for which C (1) < 1 and ANi < 1. 

An alternative manner of deriving this solution that is 
perhaps more physically intuitive, is to begin by assum- 
ing a state of the form of (|24|) and showing that the value 
of 7 that minimizes the energy is indeed the value given 
in d22). We begin by recalling the form of H 2 , exhibited 
in (H). Assuming that |£| <C 1, the T 2 term in H 2 can 
safely be neglected. If we introduce the operators 



G 



Hi = Ni — — , and 
N 

ni = N 2 - — , 

then we find that 

(#2) = £<w + (£12 + gTi(N - 1)) <j4a 2 + 4a x J> 

where Eduai is the energy of the dual condensate. Since 
the magnitude of the cross-interaction term involving T± 
only depends on the many-body state through expec- 
tation values of bilinear operators, this term together 
with the £12 term can be considered as an effective 
single-particle energy. To first order in 7 and to lead- 
ing order in 1/N, the perturbed dual condensate has 

(^a\a2 + a\a\^ = 2^N and consequently it has an ef- 
fective single-particle energy benefit over the dual con- 
densate of 2 r yN(— €12 — gNT\) (this quantity is positive 
at the barrier strengths of interest, since £12 < and typ- 
ically 7\ < 0). On the other hand, (ti\) = (h\) = 2j 2 
for such a state, corresponding to a self-interaction en- 
ergy cost of 2<y 2 gTo. Thus, the effective single-particle 
energy of the perturbed dual condensate decreases lin- 
early with 7, while the self-interaction energy increases 
quadratically with this parameter. The minimum occurs 
precisely when 7 has the value N(— 1\2 — gNTi)/2gTo, 
which approximates the value in ( J22| ) for N 1 . 

For typical double- well potentials, where Ng\Ti\ <C 
I e x 2 1 , 7 is well approximated by N{—E\2) /2gT$ to lead- 
ing order in 1/N. In this case, the transition from a com- 
pletely fragmented condensate to one that shows some 
coherence across the barrier occurs when the number of 
particles times the ratio of the inversion frequency to the 
self-interaction energy becomes non-negligible. Thus, for 
a given shape of the trap and a fixed barrier strength, 
the degree of fragmentation decreases as the number of 
particles is increased but increases as the strength of the 
interaction is increased. 



C. Numerical solutions 

The form of the double-well potential in the MIT con- 
densate interference experiment is well modelled by a 
term that is harmonic along all three Cartesian axes, with 
frequencies uj x , lu v and uj z respectively, to which is added 
a Gaussian barrier of width 5 and strength a centered at 
x = 0, 

U(r) = miulx 2 +iuy +w 2 z z 2 ) + -?L- e -& . (28) 

V Zird 

The parameter values appropriate for Ref. jlj are ui x = 
2n x 19 Hz, uj y — lo z — 2tt x 250 Hz, and 8 — 6 /im. 
Within such a trap, the criterion (^) for the applicability 



of our perturbative approach becomes TV <§; 100, which 
is much smaller than the number of condensate atoms 
in their experiment. We consider instead a larger trap, 
specifically, one which is isotropic with the trapping fre- 
quency of the axis of weakest confinement in the MIT 
trap, uj x = Lu y = lu z = 2tt x 19 Hz. In this case our 
perturbative approach is good for up to approximately 
N = 100 particles, and this is the example we consider. 
The scattering length of 23 Na is taken to be a sc = 3 nm 

ill- 

In Fig. 1 we plot the profile along the x-axis of the 
external potential U(r) and the wavefunctions <\>\ (r) and 
cj)2 (r), together with the coefficients Cjv t for several val- 
ues of the barrier strength. Fig. 2 displays the degree 
of first order spatial coherence, and the variance, 

AN±, in the number of particles occupying the localized 
state 4>i as a function of the barrier strength. The generic 
features of these results persist for different choices of pa- 
rameters in (^8|) as well as for different choices of the form 
of U(r). 

Also displayed in Fig. 2 are the values for and 
AA^! given by the continuum approximation, as specified 
by eqs. (@) and (|l|), and given by the two-coefficient 
approximation, as specified by eqs. ( p6| ) and (|27|), for the 
same choice of external potential. For their respective 
ranges of validity, the analytic approximations are found 
to fit the numerical results extremely well. 

From these calculations arise the following picture of 
the transition between a single and a fragmented con- 
densate. Moving up from zero barrier strength, there is 
a range of barrier strengths over which CW is close to 
unity, while AA^i falls from its single condensate value of 
to a value of 1. Moving down from infinite barrier 
strength, there is a range of barrier strengths over which 
CW and AA^i arc both much less than 1. Between these 
two domains, there is a narrow range of barrier strengths 
wherein the greatest part of the transition in is made. 
The barrier strengths delimiting these domains can be es- 
timated analytically using the approximations presented 
in this section. 



IV. DISCUSSION 

A. The experimental signature of fragmentation 

Herein we consider a measurement of the first-order de- 
gree of spatial coherence. This is accomplished by a type 
of interference experiment that has been widely discussed 
in the literature § |]| ||J Q (TJ. Essentially, it con- 
stitutes a double-slit experiment for Bose condensates. 
The thought experiment runs as follows. After prepara- 
tion of the condensate, the trap potential is removed and 
the atoms fall under the force of gravity through a pair of 
slits, located symmetrically about x — 0. These slits can 
be formed by changing the shape of the trapping poten- 
tial, as long as this change is not so rapid that excitations 
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are induced, and not so slow that the system has time to 
relax to a new many-body ground state. For simplicity 
of the analysis, we also assume that the particles on the 
left and right are each given momentum translations of 
magnitude hk towards one another ]T^| . In the absence 
of such translations, the interference pattern is simply 
more complicated, and has been studied by Rorhl et at 
We make the approximation that the inter-particle 
interactions are insignificant during this expansion pe- 
riod. In this case, only the single-particle wavefunctions 
evolve, while the coefficients in the many-body state 
(|J) remain unchanged. 

Suppose the slits are centered at points r = (x, y, z) 
and r' = (— x, y, z), where x is positive and chosen 
to be sufficiently large so that |</>i(r)| <C |</> 2 (r)| and 
|0i(r')| » |^ 2 (r')| ■ I n this case, the single-particle wave- 
functions (pi (r) and </> 2 (r) evolve to wavefunctions local- 
ized entirely at just one of the slits. After the momen- 
tum translation and a period of free expansion the single- 
particle wavefunctions originating from the left and right 
of the barrier acquire complex phase factors of e lkx and 
e -ikx an( j ma g n itudes we denote by 4>\(y) and </i 2 (r) re- 
spectively. For many-body states that are symmetric 
under a reflection about x = 0, these magnitudes are 
roughly uniform and of equal magnitude in the far field 
of the double-slit, so that the many-body state in the far 
field can be approximated by the many-body state prior 
to removal of the trap, with <^>i(r) and </> 2 (r) replaced by 
e ikx an( j e -ikx respectively. 

We now imagine detectors in the far field that are as- 
sumed to remove atoms from the condensate Jl6| ], The 
probability distribution over the position, ri , of the first 
detection is given by the expectation value of normally 

ordered field operators P x (ri) = jf ^ t (r 1 )^(r 1 )^) . The 
probability distribution over the positions r!,r 2 ,...,r m 
of the first m detections is given by P m (ri, r 2 , r m ) = 
^t( ri )*t( r2 )..^t( rm )^( rm )..^( ra )^( ri )^ . 

The density of detection events that emerges in a sin- 
gle run of the double-slit experiment has the form 
p m (r) = ^X^i^( r ~ r i) where the set of posi- 
tions {r^} are obtained from the probability distribution 
P m (ri, r 2 , r m ). Of course, the finite resolution of any 
realistic detector can be accounted for by replacing the 
delta function in this expression with a suitably broad- 
ened distribution; as long as the resolution is finer than 
the distance between the fringes of the interference pat- 
tern, the difference will not be significant. 

In a single run of the double-slit experiment, both the 
single and dual condensates typically yield a distribution 
p m (r) with essentially the maximum possible fringe vis- 
ibility. This is obviously true for a single condensate, 
and has been shown to be true for a dual condensate 
in the seminal paper of Javanainen and Yoo Jl6| . Thus, 
the mere presence of such interference is not indicative 
of a non-zero degree of first-order coherence. However, 
suppose the experiment is repeated many times with the 



same initial many-body state. In this case, the spatial 
phase of the interference pattern will vary randomly from 
one run to the next if the initial state is a dual conden- 
sate, while it will remain fixed if the initial state is a single 
condensate [jl6| . We therefore expect the degree of first- 
order spatial coherence to be revealed by the variance in 
the spatial phase of the interference pattern over many 
runs, or equivalently, the fringe visibility of the average 
detection pattern over many runs. Indeed, if one averages 
the pattern of detections from an infinite number of runs 
of the double-slit experiment, all prepared initially in the 
same many-body state, and each involving m detection 
events, one obtains: 

r N m 

p m {v)= / d 3 r 1 ...d 3 r m (-V<5(r-r 4 ))P m (r 1 ,r 2 ,...,r m ) 
J rn 

N m f 
m ^ J 

l — l 

xP m (r 1 , ...,ri_i,r,r i+ i, ...,r m ) 
*t( r )#( r ) N 



The final equality follows from the fact that each element 
of the sum is simply equal to P 1 (r^) . If the fringe visibility 
of the average detection pattern, p m (r), is evaluated for 
the many-body state in the far field, it is found to be 
precisely equal to our definition of CW, the degree of 
first-order spatial coherence across the barrier, and thus 
enables a measurement of the latter. 

Another possibility for an experimental study of frag- 
mentation is a measurement of the degree of second-order 
spatial coherence; this may be accessible through non- 
resonant imaging p9[. 



B. Finite temperatures 

We close this section with a few comments on the ef- 
fect of finite temperatures on the coherence properties 
of the condensate. The excited states of the system are 
not in general well approximated by a state of the form 
of (|). Nonetheless it can be shown that at infinite bar- 
rier strength the first-order degree of spatial coherence 
across the barrier is zero for thermal equilibrium at any 
finite temperature. The proof is as follows. For an in- 
finitely strong barrier the Hamiltonian is separable into 
two terms each involving only operators pertaining to 
particles on one side of the barrier. Consequently, any 
non-degenerate energy eigenstates are tensor products of 
states describing particles on one side of the barrier only. 
For any set of energy eigenstates that are degenerate, the 
subspace of Hilbert space spanned by this set has a basis 
of such product states. Thus, one can always identify a 
basis of energy eigenstates that are product states of this 
sort. Since the density operator that represents thermal 
equilibrium, p, is a mixture of these energy eigenstates, 
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we will necessarily have Tr(p^ (r)^'(r')) = if r and r' 
are on opposite sides of the barrier, and consequently the 
degree of first-order spatial coherence across the barrier 
for such a mixed state is zero as well. In the absence of 
any barrier, as long as the temperature is small enough 
that most of the particles are in the lowest single-particle 
energy level, the first-order degree of spatial coherence for 
the thermal state should be close to unity. For such tem- 
peratures, if the barrier strength is varied from zero to 
infinity, the first-order degree of spatial coherence of the 
thermal state will vary from nearly unity to zero. Thus, 
a significant transition must still occur at such tempera- 
tures. 



Our variational approach can be extended in a 
straightforward manner to the determination of the 
many-body ground state in systems where the external 
potential has an arbitrary number, n, of minima. In such 
cases, one would simply introduce states that are arbi- 
trary superpositions of Fock states where up to n single- 
particle wavefunctions are occupied. Such an analysis 
should be of relevance to the determination of the coher- 
ence properties of Bose condensates in optical lattices. 
Moreover, at extremely high densities, where the inter- 
action energy is dominant, it may become energetically 
favorable for a condensate to begin to fragment even in 
the presence of a perfectly uniform potential. We hope 
to address these possibilities in future work. 



V. CONCLUSIONS 

A theoretical treatment of fragmentation in Bose con- 
densates must go beyond a mean-field analysis. For the 
case of repulsive inter-particle interactions and a double- 
well trapping potential, we have proposed an approach 
wherein an approximation to the many-body ground 
state is obtained by a restricted variational principle. We 
have implemented this proposal for the case of nearly 
non-interacting particles. 

The coherence properties that we have considered are 
the degrees of first-order and second-order coherence 
across the central barrier of the potential. The first of 
these quantifies the variance over many runs in the spa- 
tial phase of the fringe pattern arising from the interfer- 
ence of atoms on either side of the barrier. The second is 
essentially the density-density correlation across the bar- 
rier, and for the states we consider it is a simple function 
of the variance in the number of particles in one of the 
wells. We find that as the barrier strength is increased, 
this variance is continuously squeezed down from its value 
for a single condensate. The degree of first-order spatial 
coherence is close to unity when this variance is greater 
than one, but thereafter drops off rapidly. Above a cer- 
tain critical barrier strength we find that both quantities 
become much less than one, indicating that the conden- 
sate is essentially completely fragmented. 

We have discussed how the degree of first-order co- 
herence might be measured through interference exper- 
iments, and argued that a significant effect should be 
present even at finite temperatures. A concern, however, 
is that the ground state might be difficult to prepare if 
the relaxation time of the system is long compared to the 
lifetime of the condensate. This could arise if the only 
way for the particles to be redistributed across the bar- 
rier is by tunneling through it. However for numbers of 
particles that are not too large, this tunneling time need 
not be restrictive. For instance, in the example presented 
in Section III.C, the single particle tunneling time at the 
barrier strength where = 0.88 is approximately one 
minute, while at the barrier strength where —0.1 it 
is roughly one hour. 
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Fig.l. Numerical solutions for the single-particle wave- 
functions 4>\ (r) and 02 (r) and the coefficients Cjvi for 
N = 100 particles and double- well potentials with bar- 
ri er streng ths of a = 0,15,30,45, and 60 in units of 
yxwj/m. The remaining parameter values are speci- 
fied in the text. The dotted curve is the external poten- 
tial along the x-axis, U(x, y = 0, z = 0), in units of Huj x . 
The solid and dashed curves are tf>\(x,y = 0, z = 0) and 
02 (x, y — 0, z = 0) respectively, in arbitrary units. 

Fig. 2. (a) The degree of first-order spatial coherence 
across the barrier, and (b) the variance in the oc- 
cupation number of one of the wells, AiVi, as a func- 
tion of the barrier strength, a, in units of ^JWux^jm for 
N = 100 particles and the parameter values specified in 
the text. The solid curve is the numerical solution, the 
dashed curve is the continuum approximation, and the 
dotted curve is the two- coefficient approximation. 
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